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Abstract 

We describe an algorithm for using a quantum computer to calculate mean values of 
observables and the partition function of a quantum system. Our algorithm includes 
two sub-algorithms. The first sub-algorithm is for calculating, with polynomial ef- 
ficiency, certain diagonal matrix elements of an observable. This sub-algorithm is 
performed on a quantum computer, using quantum phase estimation and tomog- 
raphy. The second sub-algorithm is for sampling a probability distribution. This 
sub- algorithm is not polynomially efficient. It can be performed either on a classical 
or a quantum computer, but a quantum computer can perform it quadratically faster. 



1 



1 Introduction 



In this paper, we describe an algorithm for using a quantum computer to calculate 
mean values of observables and the partition function of a quantum system. Our 
algorithm includes two sub-algorithms. 

One sub-algorithm is for calculating, with polynomial efficiency, certain diago- 
nal matrix elements of an observable. This sub-algorithm is performed on a quantum 
computer, using quantum phase estimation and tomography. This sub-algorithm is 
very similar to the algorithm of Ref. [T] by Harrow et al., except that we modify it 
to accomplish a substantially different job that has nothing to do whatsoever with 
systems of linear equations. 

A second sub-algorithm is for sampling a probability distribution. This sub- 
algorithm is not polynomially efficient. It can be performed either on a classical or 
a quantum computer. However, a quantum computer can perform it quadratically 
faster that a classical computer, if one uses a quantum sampling technique based 
on Szegedy operators, like, for instance, the quantum Gibbs sampling algorithm de- 
scribed in Ref. [2]. 

We end the paper with a brief section comparing the algorithms proposed in 
this paper with quantum algorithms proposed in earlier papers for calculating the 
same things. In particular, we compare our work to Ref.j3] by Wocjan et al., Ref. jl] 
by Poulin and Wocjan, and Ref. [5] by Temme et al.. 

2 Location, Loeation , Loeation 
Notation, Notation, Notation 

In this section, we will define some notation that is used throughout this paper. 
For additional information about my notational quirks, I recommend that the reader 
consult the notation section of some of my previous papers; for example, Ref. [6]. 

We will often use the symbol Nb for the number (> 1) of qubits and N$ = 2 Nb 
for the number of states with Nb qubits. The quantum computing literature often 
uses n for Nb and N for N$, but we will avoid this notation. We prefer to use n for 
the number operator |1)(1|. 

Let Bool = {0, 1}. As usual, let Z, E, C represent the set of integers (negative 
and non-negative), real numbers, and complex numbers, respectively. We will also 
sometimes add a superscript to the symbols Z, K to indicate a subset of these sets. 
For example, we will use R- to denote the non-negative reals. For integers a, b such 
that a <b, let Z a h = {a, a + 1, ... b — 1, b}. 

We will use 0(5*) to represent the "truth function"; Q(S) equals 1 if statement 
S is true and if S is false. For example, the Kronecker delta function is defined by 
5y = 5(x,y) = Q(x = y). 

If x = xn b -i ■ ■ ■ x 2 xix , where x^ 6 Bool, then dec{x) = J2u=o 1 = x. 



2 



Conversely, x = bin(x). However, when our meaning is clear from context, we will 
omit the binQ and dec(). Hence, in some places x might stand for an element of 
Zo,n s -i, and in other places for the corresponding element of BooI Nb . 

We won't usually put a caret over a symbol to indicate that it is an operator, 
but sometimes we will. For example, we will usually use H for a Hamiltonian, but 
sometimes, for clarity, we will call it H . 

Note that in our quantum circuit diagrams, time flows from the right to the left 
of the diagram (this is the Dirac Convention). Careful: Many workers in Quantum 
Computing draw their diagrams so that time flows from left to right (the Quayle 
Convention). 

We will say a problem can be solved with polynomial efficiency, or p- 
emciently for short, if its solution can be achieved in a time polynomial in Nb- Here 
Nb is the number of bits required to encode the input for the algorithm that solves 
the problem. 

By compiling a unitary matrix, we mean decomposing it into a SEO (Se- 
quence of Elementary Operators), where by elementary operators we mean operators 
that act on only a few qubits (usually 1, 2 or 3), such as single-qubit rotations and 
CNOTs. Compilations can be either exact, or approximate (within a certain preci- 
sion) . 

We will say a unitary operator U acting on C s can be compiled with 
polynomial efficiently, or p-compiled for short, if U can be expressed, either 
approximately or exactly, as a SEO of length polynomial in Nb- When necessary, we 
specify whether a p-compilation is exact or approximate. 

Next, we explain our notation related to Discrete Fourier Transforms. 

For any x G Z 0j ff s -i, define k x by 

a) 

Let k and x be operators acting on C Ns with eigenvectors and eigenvalues given by 

x\x) = x\x) (2) 

and 

k\k = k x ) = k x \k = k x ) (3) 

for any x G Z 0i jy g _ 1 . Let the eigenvectors of x and k be related by a "Discrete Fourier 
Transformation" : 

1 N S -1 

\k = kx ) = ^=Y. elkxy \y)- w 

V S y=0 

Eq.fllJ) defines a "basis-changer" unitary operator Uft with matrix elements given by 
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(y\U FT \x) = (y\k = k x ) = -^=e^ = ^=e ixky . (5) 

v v Ms 

Assume that the eigenstates of x are orthonormal and complete: 

(y\x) = 8(y,x) (6) 

for all x,y G Z 0i jv s -i, and 

N s -l 

X)|x><x| = l. (7) 

x=0 

Then it follows that the eigenstates of fc are orthonormal and complete too, because 
(k = k y \k = k x ) = (y\U FT U FT \x) = 5(y, x) , (8) 

and 

J2\k = k x )(k = k x \=J2 U FT\x)(x\U FT = l. (9) 

X X 

Often in quantum computing we come across quantum states of the form 



, Ns-l 

\k = k) = ^=J2 elky \y) > ( 10 ) 

s y=o 



with 

k = k z + Ak, (11) 

for some z G Z 0j n s _ 2 and < A A; < •j^. If A; is neither k z nor fc 2+ i but lies somewhere 

in between, then \k) is close to but not exactly equal to an eigenstate of k. Note that 
when k = k z , 

(x\U FT \k = k z ) = (x\z) = S(x,z) . (12) 
For k = k where k is given by Eq. fllll) . this generalizes to 



N S -1 



i^oc I U pvj-i | k 



N s 



D i(k-k x )y 



11 — e i(k-k x )N s 



N s 1 — e i{k-k x ) 

1 e *(fc-*-)^F sin ((Jfc - k x )^) 
W s e i{k-k x )\ sin {{k - k x )^) 



(13) 

(14) 
(15) 
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3 Diagonal Matrix Elements 



In this section, we give an algorithm for calculating certain diagonal matrix elements 
p-efficiently, using a quantum computer. This algorithm will be used as a subroutine 
in the algorithms proposed in later sections. 

The algorithm described in the proof of Claim [1] below is very similar to the 
algorithm of Ref.[T], except that we modify it to accomplish a different job. In a nut- 
shell, Ref. p] combines two operations that were separately familiar to most quantum 
computerists long before Ref. PQ: a phase estimation (PE) operation, followed by a 
quantum tomography operation. 



(x (S) |U> 



-O 



H|_| > ^N Bj bit s 



V 



|x o > N B b lts 



-|0> 



1 bit 



tomo PE 
graphy 



Figure 1: Quantum circuit used in Claim [T] to calculate the diagonal matrix element 
given by Eq.f fT6|) . 



Claim 1 Let A be an N$ dimensional Hermitian matrix with non-negative eigenval- 
ues, V an Ns dimensional unitary matrix, f : R-° — > and Xq G BooI Nb , where 
Ns = 2 3 . Assume that f is simple (that is, that it can be calculated p-efficiently). 
Assume that we know how to p-compile V and e tAAt for any At > 0. Then we can 
calculate p-efficiently the diagonal matrix element 

IM(x ) = (x Q \V^f(A)V\x ) . (16) 

proof: 

Let Nsj = 2 Nbj - be the number of states of Nsj bits, for a set of bits labeled j. 

The diagonal matrix element Eq. (I16p can be calculated by running on a quan- 
tum computer the quantum circuit shown in FigJTJ In that figure, there are three 
sets of qubits. At the top are Nsj "probe" qubits. Below the probe qubits are the 
Nb "main" qubits which are initially in state |xo). Finally, below the main qubits is 
a single "ancilla" qubit that is used in the final tomography step. 
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In a moment, we will describe the evolution of the state vector as it courses 
down this quantum circuit. But before doing so, we need to specify some of the blocks 
in Fig{T]more precisely. 

Note that 

N Si -i 

H 9Nb >\0) 9Nb > = J2 \i) ■ ( 17 ) 

As for the T box, it represents 



r 



U: 



PE 



u 



PE 



u 



PE 



tfs,~i — |_b">C7l - 



3=0 



(19) 



where 



U PE = e 



iAAt 



(20) 



(For defmiteness, some parts of Eq. (1181) assume that Nsj = 3.) Finally, the operation 
with the "half moon" nodes represents 



c — 




t) — 




Nsi-l 


4) 


= E 




3=0 



R„ 



(21) 



R; 



(note that some parts of Eq. (1211) again assume that Nsj = 3.) where Rj is defined by 



Ri 



s j c j 



(22) 



with 
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_2vrj_ 
AtN 



Sj 



(23) 



(24) 



7 is defined so that < Cj < 1 for all j G Z 0)Ns .-i. Eq.(l2TT) is an SU(2) multiplexor 
with iV^j controls. It can be p-compiled approximately (the software application 
"Multiplexor Expander" [7] can do this). 
Let 



A\A — A x ) — A X \A — A x ) 
for x G BooI Nb . Define the basis-changer unitary operator Ua by 



(25) 



(y\U A \x) = (y\A = A x ) (26) 

for x, y G BooI Nb . 

Figfj] includes at its bottom a time axis marked with notches for times 1 thru 
4. Let \^ t ) for t G Z\a denote the quantum state at those times. Let us represent 
such states by sums over 3 rows, where the first row refers to the Nsj probe bits, the 
second row refers to the Nb main bits, and the final row refers to the single ancilla 
bit. 

One has 



l*i> 



E 

|o> 



N S _ 

N s - 

2 = 



' \A = 



A X )(A 



A x \V\x ) 



(27) 







1*2} = 


AA&tj' V^A^s-l 1 

L io> 


= E E 

j=0 x=0 


" \jY- 

\A = 

. |o> 



E 7 =o \j) 



(2? 



(29) 
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I*3> 



E E 

3=0 x=0 



N S -1 

E 

x=0 



|i = A x .)(i = ^|V|x ) 



V^" 1 U'\ V J, " 1_i £ 

Z^i'=0 U / 2^=0 N Si 

\A = A X )(A = A x \V\x ) 
L |0> 



Ar s _lWs z -l 

E E 

x=0 j'=0 



\f)5( 2 ^-,A x At) 



\A = A X )(A = A x \V\x ) 
10) 



iV s -l 

E 

x=o I |o) 



\A = A X )(A = A x \V\x ) 



(30) 



(31) 



(32) 



(33) 



and 



1*4) 



Ns-l 

E 

x=0 



\j 

\A = A*)(A = A x \V\x ) 



>/t7(^)|0> + >/1-7/(^)|1> 



(34) 



AtNs ■ 

Note that we are treating the quantity A x 2?r - as if it were an integer in the range 
Zq,n s —i- This quantity is non- negative because A x and At are non-negative by 
assumption. At can be chosen small enough so that this quantity is always smaller 
than Nsj — 1. As to whether the quantity is an integer, it need not be in general. 
When it isn't, one has to do a more careful treatment (as done in Ref.pQ). By using 
the treatment given above instead of a more careful one like the one of Ref.pQ, we are 
introducing a small error in what we shall say next. 
It follows from Eq.QSD that 



0/1 

(35) 

1 Note that (A = A X ( S )| = (x^\U A so if we were actually going to observe the A of the main 
qubits, this would require that we know and apply ll\ to them and then measure (a;( s )|. However, it 
turns out that we don't need to observe/measure the main qubits to calculate so we don't need 
to know or apply U A after all. 



8 



where j f G BooI Nb i, x& G BooI Nb , G Bool. The index s G Z^ Nsam labels N sam 
samples. We can define an empirical probability distribution P x ^ by 

-i N sam 

P*j,(*,b)=N—Y t % W 8t'\ (36) 

^ sam , 
s=l 

where x G Z 0i Ar s _i and b G Bool. Then 

P x , b _(x,b) « | (A = A x \y/7f(A)V\x ) | 2 S° b + \{A = A x y\- 1 f{A)V\x^ 5\ . (37) 
Therefore, 

aw = E ^ h ) = — E > ( 38 ) 

and 



TV 

2 = S ""' 8=1 



A(6) ~ (xol^T/^^ko)^ + (x \V*[l - 1 f{A)]V\x Q )8 1 b . (39) 
Finally, we conclude that 

(x \V^f(A)V\x ) « ^ = - J- E <T • (40) 

QED 

4 Mean Values of Observables and 
Boltzmann Partition Function 

In this section, we give algorithms for calculating, under various scenarios, the ex- 
pected value of an observable Q of a quantum system, and the Boltzmann partition 
function Z = e~ l3Er for an inverse temperature (5 G where E r are the eigen- 
values of the Hamiltonian of the quantum system. 

Consider a quantum system with density matrix p and Hamiltonian H, where 
both operators act on C^ 5 . (careful: H is also used for the single qubit Hadamard 
transformation). Let 

H\H — E x ) — E X \H — E x ) (41) 
for x G BooI Nb . Define the basis-changer unitary operator Uh by 

(y\U H \x) = (y\H = E x ) (42) 
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for x, y G BooI Nb . Furthermore, consider a Hermitian operator Q acting on C Ns 
(what we call an "observable" of the quantum system) with 



Q\Q = Q x ) = Q X \Q = Q x ) (43) 

and 

(y\U n \x) = (y\n = n x ) (44) 

for x, y G BooI Nb . 

Some important scenarios that we would like to consider are 

(a) Assume that Q x for all x G BooI Nb and Uq are known. Furthermore, assume 
that we know how to p-compile Uq and e lpAt for any At > 0. Then tr(fip) can 
be calculated as follows. 

In Claim [H set V = Uq, A = p and /(£) = £ for £ G M^°. Thus 



H(x) = (x\U n pU n \x) (45) 

can be calculated p-efficiently for any x G BooI Nb . Now a quantum sampling 
algorithm like the quantum Gibbs sampling algorithm of Ref.p] can be used 
to sample the probability distribution P%{x) = p(x) given by Eq. (1451) . If this 
yields the sample points x^ for s — 1, 2, . . . N sam , then 



tr(fip) « — — n„ w • (46) 

sam 

(b) Assume that Q x for all x G BooI Nb and Uq are known. Furthermore, assume 
that we know how to p-compile Uq and e tHAt for any At > 0. Let p = e ^ 
where (3 G M-° and Z = tr(e - ^ H ). T/ien tr(Op) can be calculated as follows. 

In Claim [H set V = U n , A = H and /(£) = e - ^ for ^ G R-°. (Assume that the 
eigenvalues of H are bounded below, as is true for any physical Hamiltonian, 
and, if necessary, add a constant to H so as to make all its eigenvalues non- 
negative.) Thus 



H{x) = (x\W n e-P H U n \x) (47) 

can be calculated p-efficiently for any x G BooI Nb . Now a quantum sampling 
algorithm like the quantum Gibbs sampling algorithm of Ref . [2] can be used to 
sample p(x) given by Eq. (jl7j) . (No need to know XLM X ) since the sampling 
algorithm only uses probability ratios.) If this yields the sample points x^ for 
s = 1, 2, . . . N sam , then an estimate of the expected value of Q is again given by 
Eq.diBl. 
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(c) Assume that function g : R-° — > R-° zs simple (that is, that it can be calculated 
p- efficiently) . Assume that we know how to p-compile e lHAt for any At > 0. 
Let p = e -^- where (3 e and Z = tr(e"^). Let Z g = ti{g(H)e~^ H } . Then 
Z g and tr{g(H)p} can be calculated as follows. 

In Claim [U set V = 1, A = H and /(£) = g^e'^ for £ e (If necessary, 

add a constant to H so as to make all its eigenvalues non-negative.) Thus 



Hix) = (x\g{H)e- pH \x) (48) 

can be calculated p-efficiently for any x G BooI Nb . Now a quantum sampling 
algorithm like the quantum Gibbs sampling algorithm of Ref. [2] can be used 
to sample p(x) given by Eq. (l4"8~|) . If this yields the sample points x^ s ' for s = 
1,2,... N sam , then an estimate Z g of Z g can be obtained as follows. 



I N S am / \ 

sam i 

s=l a 



a= J_'f«» V (50) 



s=l 

To estimate tr{g(if)p}, we can use 



« = tr{^(if)p} . (51) 

Note that if : M-° -> R instead of g : JR-° ->■ R-°, one can set g = g+ - g_ 
where g± > 0. Then Z g = Z g+ — Z g _. So we can estimate Z g by first estimating 
Z g± . Since we can estimate Z g and Zi, we can estimate tr{g(H)p} using Eq. floTj) . 



5 Comparisons with the Joneses 

In this section, we compare the algorithms of this paper with algorithms in earlier 
papers for calculating the same things ("our competition"). This paper has not 
addressed the problem of finding useful estimates and bounds for the convergence 
rate and error of its algorithms. Such analysis will hopefully be done in future papers 
by me or others more capable. Lacking this analysis, the comparisons in this section 
should be taken as incomplete and preliminary. 

Suppose V is a Hilbert space and p is a density matrix acting on C s . Any 
pure state e <g> V such that trv(|\l/)(\l/|) = p is called a purification of p. 

Given a probability distribution Px(x), where x E Zq^ s -i, it is convenient to 
define a | -JPx) state by 
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= E \fp^-)\ x ) 



N S -1 



(52) 



x=0 



Note that in the algorithms (a), (b), (c), finding the needed diagonal matrix 
elements is p-efficient. So the rate determining step of these algorithms is the quantum 
sampling part, which is not p-efficient. However, if the quantum sampling is done 
using the algorithm of Ref . [2j , then the sampling converges (to a given precision) in 
C(t/|) steps, where S is the distance between the two largest eigenvalue magnitudes 
of the underlying Markov chain. 

Note also that none of the algorithms (a), (6), (c) ever physically produces a 
Boltzmann quantum state (either p = e~^ H /Z or a purification thereof). They do 
produce \\fF2) states. 




Ref. [3] by Wocjan et al. presents an algorithm for calculating the partition 
function of a quantum system with Hamiltonian H, but it requires that the eigenvalues 
of H be known a priori (as is the case when the system is classical). The more recent 
Ref. [I] by Poulin and Wocjan does not require that the eigenvalues of H be known a 
priori. However, the algorithm of Ref. [I] physically creates a sequence of Boltzmann 
quantum states (purifications of p = e~^ n /Z), each state characterized by a different 
inverse temperature, with the sequence of inverse temperatures increasing gradually, 
according to an "annealing schedule", from towards the target 0. Algorithm (c) 
also calculates the partition function of a quantum system without assuming that the 
eigenvalues of H are known a priori. But algorithm (c) never produces physically any 
Boltzmann quantum state. Nor does it introduce error at every stage of an annealing 
schedule (all its calculations are done at the target j5). 

Ref. [5] by Temme et al. gives an algorithm that physically produces a Boltz- 
mann quantum state p = e~^ H jZ. Instead of producing this p or a purification of it, 
algorithm (6) produces a | yPx) state. If one wants to use the algorithm of Temme 
et al. to get tr(fip), one still must assume the same things as (b) above (namely, 
that fl x for all x and Uq, are known, and that we know how to p-compile Uq and 
e lHAt ). If these assumptions are satisfied, one can use the algorithm of Temme et al. 
to produce U^pUq physically, and then measure \x)(x\ on it. After obtaining N sam 
samples of x, one can then use Eq. (fl6l) to estimate tr(Qp). Thus, the Temme et al. 
algorithm assumes the same things as (b) to estimate tr(Qp). A big advantage of 
(b) over Temme et al. is that (b) converges in O(^) steps, whereas Temme et al. 

converges, just like classical algorithms, in 0{\) steps. 
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